#include <stdio.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <stdlib.h>
#define ll long long
#define maxn 100+5
#define Pi 3.1415926535
#pragma GCC optimize(2)
const double eps=1e-11;

int jacobiMethord(double a[][maxn][2],double u_pre[][maxn],int n);
void martixMul(double a[][maxn],double b[][maxn],double c[][maxn],int n);//C=AxB
double Abs(double x);
void swap(int *x,int *y);
double sgn(double x);

int main(){
    double a[maxn][maxn][2],u[maxn][maxn];
    int n;
    printf("输入矩阵阶数n:");
    scanf("%d",&n);
    printf("输入%d行%d列的实对称矩阵：\n",n,n);
    for(int i=1;i<=n;++i)
        for(int j=1;j<=n;++j){
            scanf("%lf",&a[i][j][0]);
            if(i==j)
                u[i][i]=1.0;
            else
                u[i][j]=0.0;
        }
            

    int now=jacobiMethord(a,u,n);
    printf("特征值和特征向量依次为：\n");
    for(int i=1;i<=n;++i){
        printf("%10.6lf: [",a[i][i][now]);
        for(int j=1;j<=n;++j)
            printf("%10.6lf ",u[j][i]);
        printf("]\n");
    }
    return 0;
}

inline int jacobiMethord(double a[][maxn][2],double u_pre[][maxn],int n){
    int count=0;
    int pre=0,now=1;
    while (count<=1000000){
        int p,q;
        double maxAbs=-1.0;
        for(int i=1;i<=n;++i)
            for(int j=i+1;j<=n;++j)
                if(Abs(a[i][j][pre])>maxAbs){
                    maxAbs=Abs(a[i][j][pre]);
                    p=i;
                    q=j;
                }
            
        double Sin,Cos,Tan,Tan2,c;
            if(Abs(a[p][p][pre]-a[q][q][pre])>eps){
            Tan2=(2.00*a[p][q][pre])/(a[p][p][pre]-a[q][q][pre]);
            c=1.00/Tan2;
            Tan=sgn(c)/(Abs(c)+sqrt(c*c+1.00));
            Cos=1/(sqrt(1+Tan*Tan));
            Sin=Tan*Cos;
        } else {
            double phi=Pi*sgn(a[p][q][pre])/4.00;
            Sin=sin(phi);
            Cos=cos(phi);
        }
        //printf("max,sin,cos=%10.6lf,%10.6lf,%10.6lf\n",a[p][q][pre],Sin,Cos);
        double u[maxn][maxn];
        memset(u,0,sizeof(u));
        for(int i=1;i<=n;++i)
            u[i][i]=1.00;
        u[p][p]=u[q][q]=Cos;
        u[p][q]=-Sin;
        u[q][p]=Sin;

        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                a[i][j][now]=a[i][j][pre];
        a[p][p][now]=a[p][p][pre]*Cos*Cos+a[q][q][pre]*Sin*Sin+2.0*a[p][q][pre]*Sin*Cos;
        a[q][q][now]=a[p][p][pre]*Sin*Sin+a[q][q][pre]*Cos*Cos-2.0*a[p][q][pre]*Sin*Cos;
        for(int i=1;i<=n;++i)
            if(i!=p && i!=q){
                a[p][i][now]=a[i][p][now]=a[p][i][pre]*Cos+a[q][i][pre]*Sin;
                a[q][i][now]=a[i][q][now]=-a[p][i][pre]*Sin+a[q][i][pre]*Cos;
            }
        
        a[p][q][now]=a[q][p][now]=0.0;
        
        double u_pre_temp[maxn][maxn];
        for(int i=1;i<=n;++i)
            for(int j=1;j<=n;++j)
                u_pre_temp[i][j]=u_pre[i][j];
        martixMul(u_pre_temp,u,u_pre,n);

        double sum=0.0;
        for(int i=1;i<=n;++i)
            for(int j=i+1;j<=n;++j)
                sum+=Abs(a[i][j][now]);
        if(sum<(0.000001)){
            printf("已经迭代了%d次\n",count);
            return now;
        } else {
            swap(&now,&pre);
            ++count;
        }
    }
    
    
        
}

inline void martixMul(double a[][maxn],double b[][maxn],double c[][maxn],int n){
    //C=AxB;
    for(int i=1;i<=n;++i){
        for(int j=1;j<=n;++j){
            double sum=0.0;
            for(int k=1;k<=n;++k)
                sum+=a[i][k]*b[k][j];
            c[i][j]=sum;
        }
    }
}
inline double Abs(double x){
    return x<0.0?-x:x;
}

inline void swap(int *x,int *y){
    int temp=*x;
    *x=*y;
    *y=temp;
}

inline double sgn(double x){
    return x>0.0?1.0:(x<0.0?-1:0.0);
}